Charge instabilities in strongly correlated bilayer systems 
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We investigate the charge-instabilities of the Hubbard-Holstein model with two coupled layers. In 
this system the scattering processes naturally separate into contributions which are either symmetric 
or antisymmetric combinations with respect to exchange of the layers. It turns out that the short- 
range strong correlations suppress finite wave-vector nesting instabilities for both symmetries but 
favor the occurrence of phase separation in the symmetric channel. Inclusion of a sizeable long- 
range Coulomb (LRC) interaction frustrates the q = instabilities and supports the formation 
of incommensurate charge-density waves (CDW). Upon reducing doping from half-filling and for 
small electron-phonon coupling g the CDW instability first occurs in the antisymmetric channel 
but both instability lines merge with increasing g. While LRC forces always suppress the phase 
separation instability in the symmetric channel, the CDW period in the antisymmetric sector tends 
to infinity {qc 0) for sufficiently small Coulomb interaction. This feature allows for the possibility 
of singular scattering over the whole Fermi surface. We discuss possible implications of our results 
for the bilayer high-Tc cuprates. 

PACS numbers: 71.27.-|-a,74.72.-h,74.25.Kc 



I. INTRODUCTION 

Among the variety of cuprate superconductors most 
angle-resolved photoemission spectroscopy (ARPES) ex- 
periments have been made on the bilayer compound 
Bi2Sr2CaCu208+5 (Bi2212) due to the advantage of good 
cleavage planes and its nearly perfectly two-dimensional 
electronic structured However, only recent improvements 
in the resolution of ARPES measurements allowed for 
the detection of the band splitting due to coherent c- 
axis coupling of the two layers within a unit cell .^•'^ The 
same feature has also been observed in modulation-free 
(Bi,Pb)-22124 and there is experimental evidence that 
the magnitude of the bilayer splitting is constant over a 
large range of dopingj^ The splitting obeys the expected 
symmetry of LDA computationsSii being essentially zero 
along the diagonals (0, 0) (tt, tt) and (in the normal 
state) acquires a value between 88meV^ and llOmcV'^ at 
the (tt, 0) point of the Brillouin zone (BZ). 

Below Tc, the ARPES spectra display additional fea- 
tures which can be interpreted in terms of a coupling 
of the charge carriers to a collective mode. As a conse- 
quence the dispersion along the nodal direction shows a 
break at some characteristic energy uj^r' with an increased 
effective mass for binding energies smaller than loq. It has 
been shown that an analogous anomaly is also present 
in the bonding band dispersion of Bi2212 around (7r,0) 
(M-point) which reveals as an additional peak-dip hump 
feature in the ARPES line shape»i2i However, more re- 
cent ARPES experiments on underdoped (Bi,Pb)-2213ii 
have revealed a similar mass renormalization in bonding 
and antibonding band. 

Concerning the physical origin of these mode-type fea- 
tures the various proposals include a magnetic resonance 
(see e.g. Refs. 12,13 14), the coupling to phonons^ or in- 
commensurate charge-density waves (ICDW)^^'^^ as the 
source of the associated scattering. Further on from 



neutron scattering experiments (see e.g. Ref. Il7|) it 
is known that the spin fluctuations between the layers 
in YBa2CU307_5 are antiferromagnetically correlated so 
that the corresponding exchange potential is antisym- 
metric. However, since the weight of the resonance is 
quite small, large coupling constants are required in or- 
der to reproduce the spectral features within a magnetic 
mechanism and therefore this scenario is still controver- 
sially discussed.-® 

In this paper we show that antisymmetric scattering 
in a bilayer system is not unique to a magnetic interac- 
tion but may also occur in the charge sector close to a 
ICDW instability. Our analysis can be viewed as an ex- 
tension of the theory proposed by Castellani et al. in Ref. 
HqI according to which the anomalous electronic proper- 
ties of high-Tc materials are determined by a quantum 
critical point (QCP) located near optimal doping. The 
occurence of such an instability towards ICDW forma- 
tion can be theoretically substantiated by considering the 
interplay between phase separation (PS) and the long- 
range Coulomb interaction. PS is a natural feature of 
systems where strong electronic correlations lead to a 
substantial reduction of the kinetic energy. As a con- 
sequence short-range range attractive interactions (e.g. 
of phononic or magnetic origin) may dominate and in- 
duce a charge aggregation in highly doped metallic re- 
gions and a simultaneous charge depletion in spatially 
separated regions. It was pointed out in Ref. |22 that 
long-range Coulomb forces oppose the charge separation 
suppressing long-wavelength density fluctuations. It has 
been shown^^'^" that this 'frustration' of PS may re- 
sult in a finite-momentum instability corresponding to a 
ICDW quantum critical point. Near this instability the 
associated singular scattering favors the occurrence of d- 
wave superconductivity'^'^ and can explain the anomalous 
hump- type absorption in the optical conductivity of over- 
doped cupratesi^ In addition, the strong fluctuations 
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associated with the proximity to a QCP may account 
for the dependence of the pseudogap temperature on the 
characteristic time scale of the particular experimenti^i 

Further experimental evidence for the existence of a 
QCP in the phase diagram of high-Tc cuprates has been 
growing over the last few years (see e.g. Ref. and 
references therein). The idea that the associated order 
in the underdoped regime is compatible with ICDW for- 
mation is also supported by recent scanning tunneling 
microscopy (STM) experiments^ These measurements 
have revealed the existence of a non-dispersing peak in 
the Fourier transformed local density of states in slightly 
overdoped Bi2212 and thus the presence of static charge 
order in this compound (for a more detailed discussion 
on the detection of charge order in cuprates see Ref. '2^. 

Previous investigations of the ICDW-QCP 
scenarioiSiSli have been based on the two-dimensional 
electronic structure of a single Cu02 plane. However, 
real cuprate compounds can be prepared with a variable 
number of Cu02 layers per unit cell which additionally 
are electronically coupled along the c-axis. Such a 
layered structure has profound consequences on the 
momentum dependence of the long-range Coulomb 
interaction and on the spectrum of low-energy collective 
modes. In this paper we study the simplest extension 
of the single-layer case, namely we investigate possible 
charge instabilities in a system consisting of two coupled 
layers. In Sec. II we introduce the Hubbard-Holstein 
bilayer model and outline the evaluation of the relevant 
effective interactions between quasiparticles within an 
expansion. In a bilayer system these interactions 
can be either symmetric or antisymmetric with respect 
to exchange of the layers. We show in Sec. Illa that 
in the absence of long-range interactions and similar 
to the single-layer system the strong local repulsion 
favors the occurence of a phase separation instability in 
the symmetric sector. Long-range Coulomb forces are 
introduced in Sec. Illb. It turns out that in this case the 
preferred charge-instability occurs in the antisymmetric 
channel although the symmetric instability line can be 
rather close. We finally discuss possible implications of 
our results for the high-Tc cuprates in Sec. IV. Since 
this paper extends previous calculations for single-layer 
systems the reader is recommended to study Ref. ,2^ for 
more details on the subject. 



II. FORMALISM 

A. The Model 

Starting point Hamiltonian is the Hubbard-Holstein 
model for two coupled layers: 



ij,(T,OL ij.cr 



(1) 



where fl^J^ annihilates (creates) an electron at site Ri of 

layer a = 1,2 and Ui^^a = fL,afi'y,a- Note that Eq. ((TJ 
does not contain long-range forces which will be included 
in Sec. Illb. 

The chemical potential is denoted by /xq and Uj and 
tjj are hopping amplitudes in- and between the layers, re- 
spectively. In the following we take as the Fourier trans- 
formed of the interlayer hopping 



t-L(fc) = 4 [(cos(fc^) - COs{ky)f/A] 



(2) 



motivated by LDA calculations^ and ARPES 
experiments^ for high-Tc bilayer compounds. 

The operators A[^^ describe dispersionless phonons 
(frequency wq) interacting with the electrons via a local 
(Holstein-type) coupling. Note that the electron-phonon 
coupling vanishes on the mean-field level since it is incor- 
porated only via the density fluctations®. 

Furtheron we take the limit U ^ oo which can be con- 
sidered within a standard slave-boson technique. ^^'^^ In 
order to implement the constraint of no double occupa- 
tion the original fcrmion operators are decomposed as 



'-ia,abia: ficr,a ^ ^L^io-.a- Moreovcr it is Conve- 
nient to introduce the limit of large orbital degeneracy 
N to introduce a small parameter 1 /N for a perturbative 
expansion. The new fermion and boson operators are 
related by the constraint 



icra'-iCTQ I bl^bia — N/2 



(3) 



which is implemented below by introducing an additional 
local Lagrange multiplier . Within the large N expan- 
sion the model can then be represented as a functional 
integral 



Z = y DclDc^Db'^ DbDXDADA^ e~ ^0 ^"^^ 



(4) 



E4 



yb]— hA' 



dr ' dr 



dr 



H (5) 



with 



= ^ E ^i34a,aC3'y,o,b]J)^^c,+^{-^^a+iX,»)nia,a 
+ ^ E (4,lCj.,2fo],2^«,l + H.c) + ^0 E 4a^»," 

+ ^iK,cc \ blab,,c. - yj 

i.OL ^ ^ 



(6) 



where we have rescaled the hopping t^^'^ t['^'' /N and 



the el. -ph. coupling constant g — > g/\N in order to 
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compensate for the presence of N fermionic degrees of 
freedom per site. The average number of particle per 
cell and plane is n = rii = n2 = (1 — S)N/2 and S — 
corresponds to half-filling, when one half electron per cell 
and per spin flavor is present in the system. 

The mean-field self-consistency equations are obtained 
by requiring the stationarity of the mean-field free energy 
and they determine the values of — (&i,2)^ = = 



Hamiltonian reads 



Q and of Ao = (Ai,i) = (A 



Then the mean-field 



MF — 



k(T 



+ 2NlNXo r, 



(7) 



where we have transformed to the bonding/antibonding 
representation for the fermionic operators 



rA 

J ka 



nB 

Jka 



1 

1 



Cfc<T,2) 
Cfc<T,2) 



(8) 
(9) 



and the respective energies are given by E^^^ — rQ{ek ± 
t-^{ky). £k is the bare in-plane dispersion which comprises 
nearest (t) and next-nearest neighbor [t' = jt) hopping 

Ek — ~2t[cos{kx) + cos{ky) + 27C0s(/ca;) cos{ky)] (10) 

and Nl denotes the number of sites per plane. Note that 
at this level the square of the mean-field value of the 
slave-boson field Bq, 6q = Ntq — N6/2, multiplicatively 
reduces both the inter- and intralayer hopping. Fig. 1 
shows the bonding and antibonding band for selected 
cuts through the Brillouin zone. According to our choice 
for the interlayer hopping Eq. ^ the splitting is largest 
at the (tt, 0) points and reduces to along the zone di- 
agonals. As far as Aq is concerned, this quantity rigidly 
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FIG. 1: Bonding and antibonding bands with respect to the 
Fermi level for parameters 7 = —0.2, to = 0.2, and doping 
5 = 0.15. Energies are measured in units of t. 



shifts the bare chemical potential /ig as a function of dop- 
ing and is self-consistently determined by the following 
equation 



Ao 



1 



2Nl V 

k 



2Ni 



(11) 



where f{E) is the Fermi function. 

The presence of the coupling with the phonons intro- 
duces new physical effects when one considers the fluc- 
tuations of the bosonic flelds. Since only a particular 
combination a = {A^ +^)/ (2"\//V) of the phonon flelds A 
and is coupled to the fermions, it is more natural to 
use the fleld a and to integrate out the orthogonal combi- 
nation a = {A — A^)/ {2\/N). Then the quadratic action 
for the boson field a reads 



phon 



= N 



n.i,a 



n 



-^0 t 



(12) 



where we have transformed the imaginary time into Mat- 
subara frequencies. Moreover, it is convenient to work 
in the radial gauge'^^, the phase of the field hi,a — 
\/Nri^aS^pi—i(t>) is gauged away and only the modu- 
lus field Vi^a is kept, while Ai^Q. acquires a time de- 
pendence Ai_Q, K,a + dr4>i,a- Thus One can define 
two three-component fields A'^ — {Sva, SXa, «„) where 
the time- and space-dependent components are the fluc- 
tuating part of the boson fields ri^a — ro {1 + Sri^a), 
K,a = -«Ao + SXi^a and a^^a- 

Writing the Hamiltonian of coupled fermions and 
bosons as H = Hmf + -H^bos + ffint, where Hmf is the 
above mean-field Hamiltonian, which is quadratic in the 
fermionic fields, i?bos is the purely bosonic part, also in- 
cluding the terms with the a, r and A bosons appearing 
in the action |SJ) and in ffphom Ea. H12|l . iJint contains 
the fermion-boson interaction terms. The inter- and in- 
tralayer hopping terms in the bilayer U = 00 Hubbard 
model give rise to a leading order self-energy contribution 
in the quadratic part of the bosonic Hamiltonian 

S™*™(g) = ^ E ^'^-^ [/(^fe"") + /(^'f )] (13) 

- ^Y.*^{k-q)[fiEi)-f{Ei)]. 

For the following it is more convenient to transform 
also the bosonic fields to symmetric and antisymmet- 
ric combinations respectively, i.e. A±. = {A^ ± A2) 

which we combine into a single vector as ^ = (^+ , A- ) . 

After Fourier transformation to momentum space, the 
bosonic part of the action reads 

H^os = NY,A^'{q)B^^''{q)A'{~q) 
qp.11 



4 



without explicitly indicating the frequency depen- 
dence for the sake of simplicity and fijV ~ 
r-|_, A+, a+, 7'_, A_, a_. The matrix S'^'", can be ex- 
plicitely determined from Eas.lf5|)- (|13|l which results in 
the following block diagonal structure 



B 



B+ 
B 



and B^ denote the 3x3 matrices 



(14) 



The factor 2 multiplying the boson matrix B arises from 
the fact that the bosonic fields in the radial gauge are 
real and n are just fermionic bubbles with insertion of 
quasiparticle-boson vertices 



f 



IUJ„ 



(22) 



B± = 





irl 





irl 
















Wo 



(15) 



The last ingredients of our perturbation theory are the 
vertices coupling the quasiparticles to the bosons. Simi- 
lar to the bosonic fields we combine them into two three- 
component vectors Anm = (^A+^„rn, A_^„rn^ allowing us 

to write the interaction part of the Hamiltonian in the 
form 



^ E E rii.^M^ {k, q) r,i.A' ii) . 

k,q,<j n7nfiv 



where the indices (nm) label the bonding and antibond- 
ing band respectively. It turns out that only the following 
vertices give rise to a non- vanishing coupling 



-25 



A 



AB 



El 



-2.g 



Ek+i 



BA = 



(17) 



-25 



We are now in the position to evaluate the self-energy 
corrections to the boson-propagators 



which can be obtained from Dyson's equation 

n = - D°nD . 

The zero order boson propagator is 

1 



D° = — B^ 
= 2N^= 



(18) 
(19) 

(20) 



From the structure of the vertices Eqs. (|17|l one can see 
that also n acquires a block diagonal structure 



„ , n+ 
^ n 



(23) 



with being symmetric 3x3 matrices. 

Possible charge instabilities of the system can be de- 
duced from divergencies in the corresponding correla- 
tion functions or scattering amplitudesi^ The above for- 
mal scheme allows to calculate the leading-order expres- 
sions of the scattering amplitude both in the particle-hole 
channel 

Tnm-M (fc, k'- q,cu) = -^Yl ^nm (k' , -Q) ^^"^ (?, Cj) A^, {k, q) 

(24) 

and in the particle-particle channel 

T^^^.^,,{k,k'-u)^-\Y.K 
X D^^^{k~k',ujWs, 



2 ^ V 2 



k'-k 



k + k' 
2~' 



k^k' 



(25) 



so that 



It should be noted that the boson propagators are of or- 
der 1/N (cf. Eq. (1201) while the occurrence of a bare 
fermionic bubble leads to a spin summation and is there- 
fore associated with a factor N (cf. Eq. H22I) '). From Eq. 
(I21|l it thus follows that in this 1 /N approach the quasi- 
particle scattering amplitudes are residual interactions of 
order 1/N. 

Since both the boson matrix B and the polarizability 
matrix n are block diagonal the same also holds for the 
scattering amplitudes. Moreover for the scattering of two 
quasiparticles on the Fermi surface (fc = kp] k' = k'p) one 
has only two different elements for the effective scattering 
amplitude and we find in the particle-hole channel 

Fs = ^AA:AA — ^AA;BB = ^BB;BB = ^BB-.AA 
Fa = TaB-AB — ^BA;BA — ^AB;BA — ^BA;AB 



D = [2iVB + n]- 



(21) 



with 



For later use we also report here the evaluation of the 
density-density response function which for the bilayer 
system is denoted as 

Pst.nAlM - {Tp'\q,lu)p''"'{q,^lL0)) (27) 
where = Ek^a fklq,afk',7 and the indices 



{s,t,n,m} = {A^ B} refer to bonding and antibonding 
band states, respectively. The density-density response 
is most conveniently expressed via the particle-hole scat- 
tering amplitudes which yields 



kk' ^■+^ ^' I 2 ^ 



The explicit form of the zero-order bubbles P% nmilj *^) 
is reported in Eqs. (|A2IA3p and from the block-diagonal 
structure of the scattering amplitude it follows that also 
the elements of Pst,nm{q, i^) decouple into the symmetric 
and antisymmetric sector. 

Additionally we report an approximate derivation for 
the scattering amplitudes in appendix fXl which neglects 
the k-dependence of the vertices. This approach provides 
a more direct insight in the basic physical aspects of the 
problem and we refer to it in the following where appro- 
priate. 



III. RESULTS 

A. Phase Separation 

Let us first consider the system without long-range 
forces. In this case our model is characterized by strong 
on-site correlations which enable the electron-lattice cou- 
pling to drive the system towards a phase separation in- 
stability due to the strong reduction of the kinetic energy 
of the charge carriers. In the bilayer system the occur- 
rence of a phase separation instability is signalled by a 
diverging scattering amplitude Ts in the symmetric sec- 
tor (cf. Eq. t^2Q\i ) in the limit ui — 0, q ^ 0. Since for zero 
temperature and q the numerator of the fermionic 
bubbles Eq. H22I) in the symmetric sector corresponds 
to a delta-function (/ (^^^+0) - / (£^|_o) ^ H^k - ^^)) 
the sum over k-states only picks up contributions at the 
Fermi energy and as a consequence the approximative 
scheme described in appendix ^ becomes exact. Thus 
the instability criterion for PS in the bilayer system fol- 



lows from Eqs. (jA8IA9p and reads as 

l + r^(0)[iVA(0)+iVi3(0)] -0. (29) 

Here A^a/b(0) refer to the density of states of bonding 
and antibonding band at the chemical potential and 
is defined in Eq. (|A6|) . It is interesting to observe that 
to lowest order there is no influence of t-^ onto the PS 
instability. Since S™*^*" < one could conclude from Eq. 
(|A6|I that interlayer charge fluctuations work coopera- 
tively with the electron-phonon interaction and support 
the PS instability. However, interlayer hopping simul- 
taneously enhances the kinetic energy which reflects in 
an enhancement of the Lagrange parameter Aq (cf. Eq. 
(Illfl ). In fact, the self-energy contributions cancel out the 
Ao term in the 2nd order scattering amplitude which 
reads as rl{uj ^ 0,q ^ 0) = j^[-2EF/S - g^/uj„] and 
only indirectly depends on the interlayer hopping via the 
Fermi energy. Here the term proportional to the Fermi 
energy {Ep < 0) corresponds to the residual repulsion 
between the quasiparticles on the Fermi surface. Despite 
the fact that we started from a infinite on-site repulsion 
between bare particles the large screening in the system 
gives rise to a finite scattering amplitude Fg. Therefore 
the additional el. -ph. coupling can always turn the in- 
teraction into an attractive one and eventually drive the 
system towards phase separation. 

In order to analyze in more detail the q = insta- 
bilities in the symmetric and antisymmetric sector it is 
instructive to modify the model in Eq. © by coupling 
the phonons to the full electron density rather than to 
density fluctuations. In this way an el. -ph. coupling is 
effective already at the mean-field level, where a nonzero 
mean-field value of the phonon field arises in both sym- 
metric (a^ — (a+)) and antisymmetric (a'i = (fl-)) com- 
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binations. It is then straightforward to show that (a_|_) 
induces a doping dependent correction to the chemical 
potential and Eq. (|29(l is identical to the condition of 
a stationary point in signaling a divergence in the 

compressibility n = —dS/dfi, i.e. phase separation (cf. 
appendix in Ref. l2ri) . On the other hand the instability 
in the antisymmetric sector corresponds to a second order 
phase transition where the order parameter (a_ ) starts to 
acquire a finite value. As a consequence of different dis- 
tortions in the two planes also the corresponding charge 
densities will be different. 
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0.1 0.2 0.3 0.4 

doping 8 

FIG. 2: Phase diagram electron-phonon coupling g versus 
doping 5 for g = instabilities. The Maxwell construction for 
the (symmetric) phase transition separates the stable from 
the phase separated region. Parameters: t =leV, 7 = —0.2, 
tx/t ~ 0.2. The inset sketches the symmetry-broken state 
within the symmetric (i.e. PS) and antisymmetric channel. 
Shaded areas indicate enhanced charge density. 

Fig. 2 displays the phase diagram el. -ph. coupling g 
versus doping (5 together with a sketch to elucidate the 
symmetric and antisymmetric q — Q instabilities. For the 
PS instability line the van-Hove singularities of bonding 
(BB) and antibonding band (AB) reflect as the two kinks 
at (5 w 0.1 and 8 ~ 0.3, respectively. The concentration 
range in between is characterized by an AB Fermi surface 
(FS) centered around F = (0, 0) and a BB FS centered 
around X = (tTjTt). Since evaluation of the fermionic 
bubbles in the antisymmetric sector Eq. (|22l) requires 
the summation over an area which is determined by the 
difference of BB and AB Fermi surfaces, F^ is strongly 
enhanced for 0.1 < (5 < 0.3. It is therefore within this 
concentration range where the antisymmetric 2nd order 
phase transition occurs before the phase separation insta- 
bility. However, a Maxwell construction has to be done 
in order to properly determine the coexistence region in 
the symmetric sector. The corresponding phase bound- 
ary is shown by the dashed line in Fig. 2. From this 
Maxwell construction we thus conclude that also in a bi- 
layer system which is strongly susceptible to a 2nd order 
instability in the antisymmetric sector the presence of 



strong correlations favors the transition towards a phase 
separated regime. 

In principle our previous analysis does not exclude the 
occurrence of a nesting induced phase transition before 
the PS instability line is reached. In order to demon- 
strate that the instability really takes place at wave vec- 
tor (7 = we report in Fig. 3 the static scattering ampli- 
tudes in the particle- hole channel V ^(kp ^k'p^q^uj = 0) 
obtained from Eq. (|26|l for a fixed doping 8 = 0.25 and 
tuning the el. -ph. coupling towards the instabilities. In 
the absence of el. -ph. coupling {g = 0) the interaction 
between quasiparticles is repulsive (i.e. F > 0) in both 
the symmetric and antisymmetric channel. Due to the 
strong local correlations this residual repulsion increases 
with increasing wave vector q in both channels. Upon 
switching on the el. -ph. coupling Tg/j^{kp,kp,q,LL! = 0) 
therefore becomes attractive for small q and diverges at 
some critical value gcrit- Note that results in Fig. 3 are 
evaluated for doping S — 0.25 where the instability in 
the antisymmetric sector occurs first. For this reason F^i 
becomes more negative for g — > than Fg. 




1 2 3 1 2 3 



(qx.O) (q„0) 

FIG. 3: Scattering amplitudes in the symmetric (a) and anti- 
symmetric (b) channel for a scan along {qx,0). The three 
curves correspond to el. -ph. couplings g — 0.2,0.21,0.22 
respectively. Parameters: 7 — —0.2, t^/i ~ 0.2, doping 
S = 0.25. 

Let us now turn to the dynamical properties in the 
absence of long-range interactions. Since we consider a 
strongly correlated bilayer system the zero sound mode 
consists of two branches. The acoustic one in the sym- 
metric sector disperses as uJs = vfV^ + UQ for a 2-d sys- 
tem and an effective interaction u& A second, optical 
branch exists in the antisymmetric sector and its energy 
scale is determined by the interlayer hopping uja ^ t-^. 
As a consequence the Boson propagator Eq. H21|l has 
two poles for each symmetry: one is the phonon and 
the other is the zero sound (which in the symmetric sec- 
tor becomes the 2-d plasmon when long-range Coulomb 
forces are included). The two modes would cross each 
other (at rather small q ^ luq/vf) if the phonon were de- 
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coupled from the fermions but repel when the coupling 
is switched on. Consequently one observes two spectral 
features in each channel. For the symmetric combination 
(and small q): (a) the zero sound at low energy which 
upon increasing g is pushed down and, becoming softer 
and softer, drives the PS; (b) the phonon mode at energy 
higher than ujq which is hardened since it is pushed up by 
the repulsion with the zero sound. For the antisymmetric 
channel (a) at low momenta and small energies appears 
the phonon which now upon increasing g softens towards 
the 2nd order instability and (b) at higher frequencies 
the zero sound optical mode which for larger momenta is 
rapidly shifted to higher frequencies and loses intensity. 
In Fig. 4 we show the dispersion of these excitations 
along the (1, 0) direction which can be obtained from the 
poles of r5/^(fc_F, fc^, w) in Eq. if^ . 
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FIG. 4: Phonon and zero sound dispersions along the {qx,0) 
direction in the symmetric (oj"*", solid) and antisymmetric , 
dot-dashed) channel. Parameters: 7 = —0.2, t'^/t = 0.2, 
g = 0.22, uo = 0.04,doping S = 0.25. 
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B. Inclusion of long-range interactions 

The nature of possible instabilities in the system (i.e. 
phase separation or finite-q charge instabilities) crucially 
depends on the structure of the long-range Coulomb 
(LRC) potential in the bilayer system. Although we ex- 
pect the effect of LRC forces to be most effective in the 
small momentum-transfer case, where the underlying lat- 
tice structure is less visible, we explicitely take into ac- 
count the real symmetry of the bilayer square-lattice sys- 
tem (cf. Fig. 5). 

The intra- and interlayer contributions to the Coulomb 
potential are derived in appendix IbI and give rise to the 
following interaction part of the Hamiltonian 



He ^ — V V'^p'^p''. 



(30) 



where 



y^ntra^^ = 0) 

Vc 



MQx,qy) 



2k;(1-k;) 
y^nter^^ = 0) 

2 



(31) 



.(32) 



AHq^,qy)-l 



2k(1-k) 



1 



- 1 



The Coulombic couphng constant Vc = e'^az/ {2e±a1y) 
has to range from roughly 0.5-3eV in order to have holes 
in neighboring Cu02 cells repelling each other with a 
strength of 0.1-0.6 eV. Concerning the lattice parame- 
ters typical values in case of YBCO are k = d/a^ ~ 0.36 



and aj^y/az 



0.32. 



Note that V^'^ is the potential be- 



tween electrons in a two-dimensional bilayer lattice and 
the small momentum behavior reads as 



i-1 



FIG. 5: Sketch of the bilayer structure which is used for eval- 
uation of the inter- and intralayer Coulomb interaction 



yv 



yr 



Vc 1 



yintra _ 



(33) 
(34) 



Upon transforming the Coulomb potential to the 
(anti) symmetric representation 



1 



Vs/Aiq) = [VSliTil) ± VifJiq)] 



(35) 



it thus turns out that the Coulombic contribution ap- 
proaches a constant value (~ 0.1...0.5ey depending on 
parameters) in the antisymmetric sector for small mo- 
mentum transfer whereas the q — divergence in the 
symmetric part of the interaction naturally leads to a 
suppression of phase separation as we will demonstrate 
in the next section. 
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0,35 




0,2 0,4 0,6 
doping 8 



0,2 0,4 0,6 0,f 
doping 5 



FIG. 6: ICDW instability lines in the phase diagram electron- 
phonon coupling g versus doping 5. a.) Vc/t = 0.4; b) Vc/t — 
2. Parameters: 7 = -0.4, uo/t = 0.04, t'i/t = 0.2,Vc/t = 
2,6 = 6, d/az — 0.5, a^y/az = 0.32. 



out by calculating the density-density response function 
P^^{q,uj). Using the corresponding short-range part 
P^^{q,Lo) Eq. H28(l one can perform the following re- 
summation 



,LR 



SR I -nSR xrCoul t>LR 



(36) 



where the matrix V*'-'""' is given by 



V 



Coul 



/ V'^iq) V'^iq) 

V^{q) V^q) 

V'^iq) V^iq) 

V V^{q) V^{q) 1 



(37) 



Inverting Eq. yields 



C. Analysis of CDW instabilities 



>LR 



Couii-lTjSR 



(38) 



The analysis of possible instabilities under the presence 
of long-range interactions is most conveniently carried so that the instabilities can be obtained from 



J 



DET{1 - PV*^""') = = 



l-VS{q)[Pil^^{q,Q)^ 
l-V^{q)[Pl§^^B{q,Q) 



- Pl§.BAi<l,0) 



-PilABiQ.O) 



Pi§.BBi<lM 
^PliBAiQ^O)] 



(39) 



r 



for the symmetric and antisymmetric channel respec- 
tively. First it should be noted that a diverging short- 
range density-density response no longer results in a 
q = instability under the presence of long-range in- 
teractions. This is particularly obvious in the sym- 
metric channel where the Coulomb potential behaves as 
V^{q ^ 0) ^ 00. Thus one always finds a vanish- 
ing compressibility and phase separation is now ruled 
out. However, in the antisymmetric channel a diverg- 
ing short-range response P^^jy^i^ 0,0) —00 leads 
to a finite value of the corresponding long-range polariz- 
ability P„%(g = 0) ^ -l/V^{q = 0) = -2a,/{dV,) 
according to Eqs. I|33I35|I . Thus in this case LRC inter- 
actions leave the system to some extend susceptible to 
long-wavelength antisymmetric density fluctuations. 

Concerning possible instabilities under the presence 
of long-range forces the conditions in Eq. (|39|l can be 
in principle fulfilled inside the (short-range) instability 
regions of symmetric and antisymmetric channel where 
P^^((7, 0) has a positive branch up to some finite wave 
vector. Due to the fact that P^^l^^q = 0,0) stays finite 
the antisymmetric CDW instability can even occur at ar- 
bitrarily small wave vectors depending on the strength of 
the LRC interaction. 



Fig. 6 depicts the phase diagram el. -ph. coupling ver- 
sus doping for two values of the Coulomb interaction. 
Up to doping 5 ~ 0.4. ..0.5 the antisymmetric instabil- 
ity occurs at smaller coupling g than the symmetric one 
whereas both instability lines merge for larger concen- 
trations. This behavior is best understood within the 
approximate formalism given in the appendix. The cor- 
responding RPA equations for the effective interactions 
Eqs. (|A8IA9p suggest that instabilities are favored for 
those wave-vectors qcrit and dopings where (a) the resid- 
ual quasiparticle interactions Fg and F^ Eqs. l|A6IA7p 
have a minimum and (b) the 'bare' charge-charge corre- 
lations Eq. ljA2p are enhanced. Concerning (a) Fig. 8 in 
appendix ^ shows a plot of F^ and r\ for the same val- 
ues of Vc also used in the results of Fig. 6. The minima 
in these curves are determined by the relative strength of 
Coulomb interaction and the residual repulsion due to the 
slave-bosons. Since the latter part decreases with S the 
minima of F|g and F^ shift to larger qcrit when doping is 
increased (cf. inset to Fig. 8 in appendix . Moreover, 
in the limit g — > the Coulomb interaction Eq. I|32|l ap- 
proaches a constant in the antisymmetric sector and thus 
the minimum in F^ shifts to rather low momenta when 
Vc becomes sufficiently small (see the corresponding dis- 
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cussion in appendix 0. 

The value of qcrit is not only influenced by the struc- 
ture of Fg^^ but also by the 'bare' charge-charge corre- 
lation functions as mentioned above. For sizeable next- 
nearest neighbor hopping t' those are naturally enhanced 
for scattering processes connecting the high-density sec- 
tions of the (open) Fermi surface around q — (tt, 0). 
Whereas in the symmetric sector these processes are be- 
tween particle-hole (ph)states of the same Fermi surface, 
the scattering is between different (i.e. bonding and anti- 
bonding) Fermi surfaces in the antisymmetric sector im- 
plying a smaller critical wave-vector qcrit in the latter 
case. 

Summarizing, the behavior of both F^ ^ and the bub- 
bles Eq. (|A2p indicates that the critical wave- vector qcrit. 
is increasing with doping. Additionally qcrit. is expected 
to be smaller in the antisymmetric sector (where it even- 
tually tends to zero for small Vc) than in the symmet- 
ric one. Therefore, at small 5 (corresponding to a small 
critical wave- vector qcrit.) the electron-phonon coupling 
has to overcome the large 1/ qcrit Coulomb repulsion in 
the symmetric channel. However, in the antisymmetric 
sector the Coulomb repulsion approaches a constant in 
the limit g — + so that the CDW instability is reached 
for much smaller g in this case. At larger doping qcrit. 
is shifted to larger values and consequently the critical 
couplings Qcrit in both channels approximately coincide 
due to the vanishing difference between {q) and V^(q) 
for larger wave-vectors. 

Finally, Fig. 7 displays the density-density response 
functions in both channels for (5 = 0.1 close to the re- 
spective instabilities. Near the instability line the system 
is characterized by a significant quasiparticle attraction 
within a large portion of momentum space. The orien- 
tation of the critical wave- vectors is strongly determined 
by the structure of the bubbles Eq. I|A2|I . For our choice 
of t' /t = —0.4 (as appropriate for Bi2212) those are nat- 
urally enhanced along the (tt, 0) axis of the Brillouin zone 
favoring singular scattering in the same direction for both 
channels. 



IV. DISCUSSION AND CONCLUSIONS 

A. Influence on the superconducting gap structure 
in a bilayer system 

As shown in Sec. II the singular scattering ampli- 
tudes in the particle-hole channel which occur close to 
the CDW instabilities are naturally connected with sin- 
gular attractive interactions in the particle-particle chan- 
nel. In this section we briefly discuss the consequences 
for the superconducting gap structure in a bilayer ma- 
terial based on the results derived above. When we fix 
the electron-phonon coupling constant g to some mate- 
rial specific value Fig. 6 suggests the investigation of the 
two following cases. 




FIG. 7: Density-density response function in the symmetric 
(a) and antisymmetric (b) sector close to the ICDW instabil- 
ity. Parameters: 7 = -0.4, g/t = 0.23, S = 0.1, coo/t = 0.04, 
t'i/t = 0.2,Vc A = 2,e = 6, d/a^ = 0.5, a^y/a, = 0.32. 

1. For small Vc and arbitrary g (Fig. 6a) or large 
Vc and large g (Fig. 6b) the instabilities in the 
symmetric and antisymmetric sectors occur almost 
at the same doping concentration. Therefore the 
associated critical fluctuations in both channels are 
comparable. 

2. At small g and large Vc the pairing interaction oc- 
curs dominantly in the antisymmetric channel (Fig. 
6b). 

Both cases can be treated within the two-band model 
(cf. Suhl et ali^) when we restrict on the static part 

of the effective interactions Fg^"^ in the symmetric and 
antisymmetric channel. A similar analysis for the one 
band-model has been performed in Ref. The SC gap 
for bonding and antibonding band can be introduced via 

^E[r^(/^+'.T/"l-.i)+r^(/"\,T/-T,a)](40) 

Q 

and for simplicity we neglect interband pairing « 
0. In order to derive some analytical results we further 
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assume the interaction F'^/"^ to be constant within some 
energy range 2ujc around E^^. Then the self-consistency 
equation for the bilayer system reads as 

[1 - T^NaFa][1 - T^NbFb] = {T^ fNANBFAFB (42) 

where N(^a)b denotes the DOS for the 
(anti)bonding band near Ep and Fa/ b{^a/ b) = 

Let us now consider the first case where the symmetric 
and antisymmetric instabihty Unes are cfose so that both 
and display singular behavior near some critical 
doping, i.e. F'^ w F"^ ^ pmtro^2. As a consequence Eq. 
(gg simplifies to 1 = r™*™F(0)(iVA + Nb)/2 when we 
are at the transition temperature which therefore yields 

kT, = 1.14c., exp(-l/(F™*™^^^^)). (43) 

Moreover, below we have the same gap value in bond- 
ing and antibonding band = A^. Close to the insta- 
bility line the layers appear to be decoupled with respect 
to the pairing interaction which is almost completely due 
to intralayer scattering. Thus the situation in this case is 
equivalent to the single layer model investigated in Ref. 
ia with an effective DOS ^^+^^^ . On the other hand 
in the second case where symmetric and antisymmetric 
instability lines are well separated only F'^ becomes sin- 
gular near the critical doping. Note that in addition to 
a negatively diverging in-plane response this implies a 
large repulsive inter-plane interaction as can be deduced 
from Eqs. (|A8IA9p . In this case Eq. reduces to 

1 = {T^)"^ NaNbFaFb and the transition temperature is 
obtained as 

kT, = 1.14c., exp(-l/(F^VAU^)). (44) 

Below Tc we have now Ayi ^ A^ and it turns out from 
Eqs. H40I41|) that the superconducting gap in the anti- 
bonding band is determined by the pair correlations of 
the bonding band and vice versa. 

In optimally doped bilayer cuprates the van-Hove sin- 
gularity of the AB is quite close to the Fermi level. There- 
fore one should expect that around 5opt the pair correla- 
tions in the AB exceed those of the BB, i.e. {f^\f^\ |) > 

ifkjf^ll)- This implies that around optimal doping 
Ab > A. A when the scattering is dominantly antisym- 
metric but Ab = A^ when the pairing interaction takes 
place in the symmetric channel. Up to now the energy 
gaps below Tj, for both antibonding and bonding band 
have been examined in detail only for an underdoped 
modulation-free Pb-Bi2212 sample in Ref. |3 It turns 
out that both gaps are identical and deviate significantly 
from d-wave symmetry around the nodal direction. This 
may be associated with a normal state contribution to 
the gap (pseudogap) in this sample since a preliminary 
analysis revealed a similar anisotropy above T,. Within 
our analysis two identical energy gaps would imply an un- 
derlying interaction where symmetric and antisymmetric 



components are of similar strength. However, pair cor- 
relations can only be significantly influenced by the van 
Hove singularity when it is separated from Ep within an 
energy scale of T,. This is probably not the case for this 
particular underdoped sample. In this regard it would be 
interesting to repeat the same analysis for an optimally 
doped sample . 

B. Influence on the normal state resistivity 

Our last point concerns the result that in the anti- 
symmetric sector close to the instability singular fluctu- 
ations with g « are possible. In fact this feature is not 
restricted to bilayer materials but should also occur in 
single layer compounds when one includes the Coulomb 
interaction between individual layers. Denoting the in- 
plane momenta with q^^ and the perpendicular momen- 
tum with q± it is known-— that the Coulomb potential 
diverges ~ l/qj^ for q± = 0. For finite q± the potential 

^c(9|| — 0) remains finite with the smallest repulsion 
for q± = TT. Within the same model investigated in this 
paper but extended to a real three-dimensional layered 
structure the possibility arises for singular in-plane fiuc- 
tuations g|| = for q± = n. This on the other hand could 
have important consequences for temperature dependent 
transport properties such as the electrical conductivity. 
Based on a Boltzmann-equation approach Hlubina and 
Rice2& have evaluated the resistivity for models which 
are characterized by strong (critical) scattering between 
selected points on the Fermi surface. The associated elec- 
tron lifetime and resistivity displays an anomalous tem- 
perature dependence which, however, is short-circuited 
by the remainding electrons on the rest of the Fermi sur- 
face ('cold regions') which scattering rate is of the stan- 
dard Fermi liquid form l/r ~ T^. As a consequence it 
was found that the resistivity has the standard Fermi 
liquid form p ^ up to some energy scale which is 
determined by the distance to the critical point. In a 
model with singular scattering at low momenta as dis- 
cussed above all points on the Fermi surface would cor- 
respond to 'hot spots', therefore the short-circuit prob- 
lem would be prevented and the critical scattering would 
determine the temperature dependence of the resistivity 
down to T = 0. 



C. Conclusion 

We have investigated the possible charge instabilities of 
a bilayer Hubbard-Holstein model. In particular we have 
focused on the question wether these instabilities prefer- 
ably occur in the symmetric or antisymmetric channel 
with respect to the exchange of the layers. 

In the absence of long-range Coulomb interactions and 
similar to the single-layer case^S our calculations support 
the existence of phase separation arising from the at- 
tractive electron-phonon interaction. However, both the 
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symmetric and antisymmetric instability lines are rather 
close (cf. Fig. 2) and we find that also the antisymmet- 
ric symmetry-breaking (which corresponds to a second 
order phase transition) occurs at wave-vector q — (cf. 
Fig. 3b) due to the strong on-site correlations. Hence 
in the bilayer model phase separation is solely supported 
due to the Maxwell construction which only applies in 
the symmetric sector where the phase transition is first 
order. 

Inclusion of long-range forces spoils phase separation 
but finite-momentum instabilities still take place in the 
symmetric sector of the charge-charge correlations. In 
the antisymmetric sector the critical wave-vector cru- 
cially depends on the strength of the long-range Coulomb 
interaction Vc and for sufficiently small Vc and low dop- 
ing can be still around 9 « 0. Moreover, since both types 
of instabilities now correspond to a 2nd order phase tran- 
sition the Maxwell construction does not apply and one 
finds that the antisymmetric instability is now favored 
especially at low doping. 

We have discussed the above findings in the context 
of high-Tc superconductors. Recent progress in the res- 
olution of ARPES experiments has made it possible to 
separably detect the superconducting gaps in bonding 
and antibonding band respectively. We have argued that 
from the relative sizes of the gaps one can in principle de- 
duce the symmetry of the underlying interaction. Within 
the Hubbard-Holstein model we find two possibilities, 
depending on the strength of electron-phonon coupling 
g and the long-range Coulomb interaction Vc- Sizeable 
el.-ph. couphng g « 0.3t (Fig. 6a, b) imphes that in 
the quantum critical region (which in the g — S phase 
diagram is at rather large doping) symmetric and an- 
tisymmetric fluctuations are comparable. On the other 
hand the scattering is dominantly antisymmetric in case 
of small g and large Vc (Fig. 6b) when one approaches 
the instability. Moreover this antisymmetric transition 
now occurs at concentrations 6 w 0.1. ..0.2 (cf. Fig. 6b) 
which covers the range where Tc is largest in the high-Tc 
cuprates and therefore is more compatible with the quan- 
tum critical point scenario. Having in mind that antifer- 
romagnetic fluctuations in the bilayer high-Tc cuprates 
are also antisymmetric with respect to exchangge of the 
layersii both charge and spin fluctuations may easily co- 
exist and determine cooperatively the unusual properties 
of the cuprates. Unfortunately our leading order analy- 
sis in 1 /N of the U ^ 00 Hubbard-Holstein model only 
captures the charge instabilities of the model whereas an- 
tiferromagnetic correlations would appear at higher order 
in 1/A^. Therefore a long way is still to be followed in 
order to formalize the interplay between charge and spin 
degrees of freedom and to answer the question how charge 
instabilities are mirrored in the spin criticality. This in- 
triguing but difflcult issue is deflnitely beyond the scope 
of the present paper but should be definitely investigated 
in future work. 
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APPENDIX A: APPROXIMATE EVALUATION 
OF THE INSTABILITIES 

For a first analysis of the instabilities where only the 
uj = behavior of the bubbles is relevant, it is convenient 
to simplify the formalism of Sec. II by approximating the 
vertices in the following way: 



(Al) 



i.e. restricting the quasiparticle momenta to kp. As a 
consequence the structure of n Eq. H22() simplifies to 

n'-(<Z, ^„) = Ia^, (kp, q) Al {kF, ~q) PO 
where Pg^.ts are now the usual fcrmionic bubbles 



f El 



k ^"T 



El 



lOJri 



(A2) 



with s = B and t — B respectively. The matrix 
representation of Pst-^ts thus acquires a block diagonal 
structure 
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(A3) 



Eq. (|19|l can now be rewritten as a Dyson equation for 
the scattering amplitudes 



pOpOp 



and the second order scattering amplitude F^j „^ 
block diagonal 
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(A4) 
is also 

(A5) 



with only two different elements 
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The small g behavior is dominated by the LRC interac- 
tion which diverges as 1/g in the symmetric channel but 
approaches a constant in case of the antisymmetric po- 
tential. For large wave- vectors both curves merge since in 
this regime the interaction is determined by the residual 
repulsion mediated by the slave-bosons. 



Within this framework long-range Coulomb interac- 
tions can be easily incorporated by adding their sym- 
metric and antisymmetric combinations Eqs. (|35|) to the 
scattering amplitudes Eqs. (|A6IA7|I . Fig. ^ displays T^g 



and r° along the (1,0) direction of the Brillouin zone 



Finally, due to the block diagonal structure of both 
the scattering amplitude and Pfj.n^ the RPA problem 
decouples into two 2x2 matrix equations. The RPA 
scattering amplitudes for intra- and interlayer scattering 
are given by 



J 



SV-' AA:AA 
-L S 



i-ro(po 



Sy-" AA-AA 



BB-.BBi 



BB-.BBI 



^ A 



1 pO /- pO -L po 



AV^ AB;BA 
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(A8) 
(A9) 



BA:AB/ 



For w = the fermionic bubbles only display a weak 
momentum dependence up to the Fermi wave- vector kp. 
Therefore the instability vectors Qc as arising from Eqs. 
(jA8IA9p are approximately determined by the minimum 
Ailmin) of the Scattering amplitudes Eqs. HA6IA7|I . 
The insets to Fig. 8 display the doping dependence of 
kmini for two valucs of the Coulomb potential. Within 
our 2-dimensional model the minimum g„iin in T'^ al- 
ways occurs at finite (but arbitrarily small) momenta 
since Vq ^ C — Aq and the residual repulsion of the 
slave-bosons behaves as ~ q^. Note that a complete 3- 
dimensional treatment would yield V(j ^ C—Aq^ so that 
in this case a true q = Q instability could be realized. 



APPENDIX B: DERIVATION OF THE 
COULOMB POTENTIAL IN A BILAYER 
SYSTEM 

In order to derive an explicit expression for the 
Coulomb potential in the spirit of the point-charge ap- 
proximation we start from the discretized form of the 
Laplace equation. Moreover, since we assume that our 
two-dimensional model represents planes of a truly three- 
dimensional lattice we also include a third spacial dimen- 
sion. For clarity we restore the explicit dependence of the 
in-plane lattice spacing axy, which in Sec. II was set to 
unity in the square two-dimensional lattice. In the third 
space direction, instead, we assume the unit cell to have a 
lattice spacing ■ In addition each unit cell contains two 
layers separated by spacing d. Note that in the present 
case (with non-equidistant sampling points along the z- 
direction) the second derivative of a function / at point 
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FIG. 8: Scattering amplitudes along the {qx,qy ~ 0) 

direction of the Brillouin zone for Vc ~ OAt (a) and Vc — 2t 
(b). The electron-phonon interaction leads to a constant shift 
of the curves which are plotted for g — 0. The insets show the 
doping dependent length of the wave- vector corresponding to 
the minimum of Fg and F^i respectively. Parameters: t = 
ley, 7 = -0.4, S = 0.2, t^i/t = 0.2,VcA = 2,? = 6, d/a^ = 
0.5, axy/a^ = 0.32. 



Xi can be represented as 

+ hif{xi-l) - {hi + hi^i)f{xi)] 

where hi is the distance between points Xi and x.^+i. 

Due to the presence of two layers in the unit cell one 
obtains the following two coupled Laplace equations for 



13 



the Coulomb potential 



-eSiRl - 



r,= 



n=x,y 



2e± 



[d^{R}+, - R 



da^^a^ - d) 
(a, - d)cj){Rl - R]) - a,0(i?2 _ i?2)] 



ri—x,y 

where i?f denote lattice sites on plane a and e_L and e|| 
are the high-frequency dielectric constants perpendicu- 
larly and along the planes respectively. The correspond- 
ing Fourier transformed equations read as 

-e^4 = A{q^,qy)4>'^^ + [Kexp(ig^a^) + I - k] (j)^^ 



2e± 



^ A{qx,qy)(l)^ + [Ke-Kp{-iq, 



1 



with e = e||/ej^, k = d/az and the in-plane momentum 
dependence is contained in 



A{qx,qy) 



_ k(1 — k) 



[cos{axyqx) + cos{axyqy) - 2]-l 



We thus obtain for the in- and intra-plane LRC poten- 
tial in three-dimensional momentum space 



lintra 



A{qx,qy) 



4k(1-k) sm 2 J 

ea^ K exp(— zgajCz) + 1 — K 



Since we are interested in the effects of the Coulomb po- 
tential on the two-layer system, we now transform from 
qz to real space for the z = Q unit cell obtaining 



A{qx,qy) 



4e^ 



4e± 



2k(1-k) ^ 



2k;(1-k) 



- 1 



In the limit k — 1/2 and Cz ^ 2d one recovers the re- 
sult of the single-layer calculation of F. Becca et ali^ in 
which case 0™*'^'' denotes the potential between succes- 
sive layers. 
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